//
//  TestNSESolver.cpp
//  NSESolver
//
//  Created by Patrick Zulian on 3/22/11.
//



#include "TestNSESolver.h"
#include "Types.h"

#include "NavierStokesSolver.h"
#include "CFDProblem.h"
#include "BoundaryHandler.h"
#include "InitialValue.h"
#include "Obstacle.h"

using namespace nsesolver;
using namespace nsesolver::mpi;

CPPUNIT_TEST_SUITE_REGISTRATION( TestNSESolver );

TestNSESolver::TestNSESolver()
{
	//TODO
}

TestNSESolver::~TestNSESolver()
{
	//TODO
}


void 
TestNSESolver::testRHSIntegrity()
{
    const int imax = 10, jmax = 10;
    const double a = 1, b = 1, alpha = 1.0, 
    gx = 0, gy = 0, tau = 0.01, ReInv = 1.0/100;
    
    NavierStokesSolver solver;
    
    Vector2d swCorner, neCorner;
    swCorner[0] = 0; swCorner[1] = 0;
    neCorner[0] = a; neCorner[1] = b;
    
    CartesianTopology cart(2);
    CFDProblem problem(imax, jmax, alpha, swCorner, neCorner, BoundaryHandler(cart, Obstacle(imax, jmax)), InitialValue(), cart);
    Matrixd & u(problem.U());
    Matrixd & v(problem.V());
    
    u.allSet(0);
    v.allSet(0);
    
    Matrixd F(u), G(v); //Copy size and values
    
    
    
    problem.setBoundary();
    solver.computeFandG(F, G, problem, tau, ReInv, gx, gy);
    
    u(Range(1, imax), Range(1, jmax+1)) += tau*gx;
    v(Range(1, imax+1), Range(1, jmax)) += tau*gy;
    
    
    CPPUNIT_ASSERT_EQUAL( u , F );
    CPPUNIT_ASSERT_EQUAL( v , G );
    
    
    
    u.allSet(1);
    v.allSet(1);
    
    problem.setBoundary();
    
    solver.computeFandG(F, G, problem, tau, ReInv, gx, gy);
    
    
    u(Range(1, imax), Range(1, jmax+1)) += tau*gx;
    v(Range(1, imax+1), Range(1, jmax)) += tau*gy;
    
    
    Matrixd uTest(u(Range(2, imax-1), Range(2, jmax-1))), //Only internal nodes unconnected with the Boundary
    FTest(F(Range(2, imax-1), Range(2, jmax-1)));
    
    Matrixd vTest(v(Range(2, imax-1), Range(2, jmax-1))), 
    GTest(G(Range(2, imax-1), Range(2, jmax-1)));
    
    
    CPPUNIT_ASSERT_EQUAL( uTest , FTest );
    CPPUNIT_ASSERT_EQUAL( vTest , GTest );
    
    
}


void 
TestNSESolver::testLaplaceOperatorIntegrity()
{
    
    const int imax = 10, jmax = 10;
    const double a = 1, b = 1, alpha = 1.0;//, omega = 1.7;
    
    Vector2d swCorner, neCorner;
    swCorner[0] = 0; swCorner[1] = 0;
    neCorner[0] = a; neCorner[1] = b;
    
    CartesianTopology cart(2);
    CFDProblem problem(imax, jmax, alpha, swCorner, neCorner, BoundaryHandler(cart, Obstacle(imax, jmax)), InitialValue(), cart);
    
    problem.P().allSet(1);
    
    for(int i = 1; i <= problem.IMax(); ++i) {
        for(int j = 1; j <= problem.JMax(); ++j) {
            CPPUNIT_ASSERT_DOUBLES_EQUAL( 0, problem.d2pdx2(i, j) + problem.d2pdy2(i,j), Numerics::eps< double >() );
        }
    }
    
    problem.P().allSet(0);
    for(int i = 1; i <= problem.IMax(); ++i) {
        for(int j = 1; j <= problem.JMax(); ++j) {
            CPPUNIT_ASSERT_DOUBLES_EQUAL( 0, problem.d2pdx2(i, j) + problem.d2pdy2(i,j), Numerics::eps< double >() );
        }
    }
    
    
    problem.P().allSet(utils::Rand()*1000);
    for(int i = 1; i <= problem.IMax(); ++i) {
        for(int j = 1; j <= problem.JMax(); ++j) {
            CPPUNIT_ASSERT_DOUBLES_EQUAL( 0, problem.d2pdx2(i, j) + problem.d2pdy2(i,j), Numerics::eps< double >() );
        }
    }
    
}

void
TestNSESolver::testObstacle()
{
//    const double gx = 0, gy = 0;
//    const double omega = 1.7;
//    const double t0 = 0, tEnd = 1, Re = 1000;
    const double alpha = 1;
    
    
    Vector2d swCorner, neCorner;
    swCorner[0] = 0; swCorner[1] = 0;
    neCorner[0] = 1; neCorner[1] = 1;
    
    
    
    
    
    const int imax = 60; 
    const int jmax = 60;     
    CartesianTopology cart(2);
    
    //    Reference< CFDProblem > problem(CFDProblemFactory::createProblem("/Users/patrick/Desktop/Master/2nd_semester/Computational_Fluid_Dynamics/nsesolver/obstacles/_CFD_Course__karman.e",  cart, imax, jmax));
    //    
    
    Obstacle obs( imax, jmax );
    obs.addCell(10, 10);
    obs.addCell(10, 11);
    obs.addCell(10, 12);
    obs.addCell(11, 10);
    obs.addCell(11, 11);
    obs.addCell(11, 12);
    obs.assignFlags();
    
    BoundaryHandler bc(cart, obs);
    //bc.setUInflowOutflow(inflow);
    //bc.setVBoundaryType( West, MovingBelt );
    //bc.setUBoundaryType( North, MovingBelt );
    CFDProblem problem( imax, jmax, alpha, swCorner, neCorner, bc, InitialValue( InitToOnes, InitToOnes ), cart );
    
}

void
TestNSESolver::testObstacleFlags()
{
    using namespace nsesolver;
    
    Obstacle obs( 150, 150 );
    
    int x=80;
    int y0=10;
    int y=y0;
    
    obs.addCell(x, ++y);
    obs.addCell(x, ++y);
    
    ++x;
    y=y0;
    
    obs.addCell(x, ++y);
    obs.addCell(x, ++y);
    obs.addCell(x, ++y);
    
    
    ++x;
    y=y0;
    
    obs.addCell(x, ++y);
    obs.addCell(x, ++y);
    obs.addCell(x, ++y);
    obs.addCell(x, ++y);
    
    
    ++x;
    y=y0;
    
    obs.addCell(x, ++y);
    obs.addCell(x, ++y);
    obs.addCell(x, ++y);
    obs.addCell(x, ++y);
    
    obs.assignFlags();
    
    for (Obstacle::ConstObsCellIter it=obs.obsCellBegin(); it!=obs.obsCellEnd(); ++it) {
        const int i = it->i, j = it->j;
        const ObstacleCellType cellType = it->type;
        
        switch (i) {
            case 80:
            {
                //    obs.addCell(x, ++y, OSouthWest);
                //    obs.addCell(x, ++y, ONorthWest);
                switch (j) {
                    case 11:
                        CPPUNIT_ASSERT_EQUAL(cellType, OSouthWest);
                        break;
                    case 12:
                        CPPUNIT_ASSERT_EQUAL(cellType,ONorthWest);
                        break;
                }
                break;   
            }
                
            case 81:
            {
                //    obs.addCell(x, ++y, OSouth);
                //    obs.addCell(x, ++y, ONull);
                //    obs.addCell(x, ++y, ONorthWest);
                switch (j) {
                    case 11:
                        CPPUNIT_ASSERT_EQUAL(cellType, OSouth);
                        break;
                    case 12:
                        CPPUNIT_ASSERT_EQUAL(cellType, ONull);
                        break;
                    case 13:
                        CPPUNIT_ASSERT_EQUAL(cellType, ONorthWest);
                        break;
                }
                break;   
            }
                
            case 82:
            {
                //    obs.addCell(x, ++y, OSouth);
                //    obs.addCell(x, ++y, ONull);
                //    obs.addCell(x, ++y, ONull);
                //    obs.addCell(x, ++y, ONorthWest);
                switch (j) {
                    case 11:
                        CPPUNIT_ASSERT_EQUAL(cellType,OSouth);
                        break;
                    case 12:
                    case 13:
                        CPPUNIT_ASSERT_EQUAL(cellType,ONull);
                        break;
                    case 14:
                        CPPUNIT_ASSERT_EQUAL(cellType,ONorthWest);
                        break;
                }
                break;   
            }
                
            case 83:
            {
                //    obs.addCell(x, ++y, OSouthEast);
                //    obs.addCell(x, ++y, OEast);
                //    obs.addCell(x, ++y, OEast);
                //    obs.addCell(x, ++y, ONorthEast);
                switch (j) {
                    case 11:
                        CPPUNIT_ASSERT_EQUAL(cellType,OSouthEast);
                        break;
                    case 12:
                    case 13:
                        CPPUNIT_ASSERT_EQUAL(cellType,OEast);
                        break;
                    case 14:
                        CPPUNIT_ASSERT_EQUAL(cellType,ONorthEast);
                        break;
                }
                break;   
            }
                
        }   
    }
}

void 
TestNSESolver::testBoundaryValuesIntegrity()
{
    const int imax = 10, jmax = 10;
    const double a = 1, b = 1, alpha = 1.0, omega = 1.7,
    gx = 3.0, gy = 1.5, tau = 0.01, ReInv = 1.0/100;
    
    NavierStokesSolver solver;
    
    Vector2d swCorner, neCorner;
    swCorner[0] = 0; swCorner[1] = 0;
    neCorner[0] = a; neCorner[1] = b;
    
    CartesianTopology cart(2);
    CFDProblem problem(imax, jmax, alpha, swCorner, neCorner, BoundaryHandler(cart, Obstacle( imax, jmax )), InitialValue(), cart);
    Matrixd & u(problem.U());
    Matrixd & v(problem.V());
    Matrixd & p(problem.P());
    
    Matrixd rhs(imax+1, jmax+1), F(u), G(v);;
    
    u.allSet(1);
    v.allSet(1);
    p.allSet(1);
    
    
    
    
    //Just move the data around a bit
    {
        problem.setBoundary();
        
        //Compute F^n, G^n by means of (26) and (27) and Set the boundaries on F, G by means of (35)
        solver.computeFandG(F, G, problem, tau, ReInv, gx, gy );
        
        //Compute rhs for pressure equation by means of the right hans side of equation (28)
        solver.computeRHS(rhs, problem, F, G, tau);
        
        //Set p^n as the initial value for the SOR scheme
        //Done internally in SORStep(...)
        solver.SORSolve(problem, rhs, omega);
        
        //compute u_n+1 v_n+1 by (24) (25)
        solver.updateVelocities(problem, G, F, tau);
    }
    
    
    //?? do we have to reset them here ??
    
    static const double eps = 10e-9; //anyhow works with this approx order
    
    //for p
    for (int j = 1; j <= jmax; j++) {
        
        CPPUNIT_ASSERT_DOUBLES_EQUAL(p(0, j), p(1, j), eps );
        CPPUNIT_ASSERT_DOUBLES_EQUAL(p(imax+1, j), p(imax, j), eps );
    }
    
    for (int i = 1; i <= imax; i++) {
        CPPUNIT_ASSERT_DOUBLES_EQUAL( p(i, 0), p(i, 1), eps ); 
        CPPUNIT_ASSERT_DOUBLES_EQUAL( p(i, jmax+1), p(i, jmax), eps );
    }
    
    
    //for u and v
    for (int j = 1; j <= jmax; j++) {
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, u(0, j), eps );
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, u(imax, j), eps );
    }
    
    for (int i = 1; i <= imax; i++) {
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, v(i, 0), eps ); 
        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, v(i, jmax), eps );
    }
    
    //for F and G
    for (int j = 1; j <= jmax; j++) {
        CPPUNIT_ASSERT_DOUBLES_EQUAL( u(0, j), F(0, j), eps );
        CPPUNIT_ASSERT_DOUBLES_EQUAL( u(imax, j), F(imax, j), eps );
    }
    
    for (int i = 1; i <= imax; i++) {
        CPPUNIT_ASSERT_DOUBLES_EQUAL( v(i, 0), G(i, 0), eps ); 
        CPPUNIT_ASSERT_DOUBLES_EQUAL( v(i, jmax), G(i, jmax), eps );
    }
    
}
